A multiscale theory for spreading and migration of adhesion-reinforced mesenchymal cells

We present a chemomechanical whole-cell theory for the spreading and migration dynamics of mesenchymal cells that can actively reinforce their adhesion to an underlying viscoelastic substrate as a function of its stiffness. Our multiscale model couples the adhesion reinforcement effect at the subcellular scale with the nonlinear mechanics of the nucleus–cytoskeletal network complex at the cellular scale to explain the concurrent monotonic area–stiffness and non-monotonic speed–stiffness relationships observed in experiments: we consider that large cell spreading on stiff substrates flattens the nucleus, increasing the viscous drag force on it. The resulting force balance dictates a reduction in the migration speed on stiff substrates. We also reproduce the experimental influence of the substrate viscosity on the cell spreading area and migration speed by elucidating how the viscosity may either maintain adhesion reinforcement or prevent it depending on the substrate stiffness. Additionally, our model captures the experimental directed migration behaviour of the adhesion-reinforced cells along a stiffness gradient, known as durotaxis, as well as up or down a viscosity gradient (viscotaxis or anti-viscotaxis), the cell moving towards an optimal viscosity in either case. Overall, our theory explains the intertwined mechanics of the cell spreading, migration speed and direction in the presence of the molecular adhesion reinforcement mechanism.


Introduction
Cell motion plays a crucial role in various biological processes ranging from tissue formation to metastasis.The mechanical properties of the extracellular matrix (ECM) are critical environmental cues that orchestrate the spreading and migration of adherent (mesenchymal) cells [1][2][3][4][5][6][7][8].Focal adhesion (FA) sites play a key role in cell mechanosensitivity by mediating the chemical and mechanical interactions between the cell and ECM.The FA dynamics can directly impact the chemically driven localized protrusions and contractions across the cell, governed by the tightly coupled spatio-temporal distributions of the proteins, such as the Rho GTPases, ROCK or CDC42, which influence the migration dynamics [9][10][11][12][13][14].Thus, a rigorous mathematical description of cell migration must include accurate modelling of the adhesion forces at different FA sites and their coupling with the cytoskeletal dynamics driven by the intracellular signalling pathways.
Although some cells such as the U-251MG glioblastoma cells exhibit a biphasic dependence of the traction force on the rigidity of the underlying matrix [15,16], which can be explained by the classical motor-clutch theory [17,18], many other systems, such as endothelial cells and fibroblasts, display a monotonic rigidity-force relationship [19,20].This behaviour is attributed to active adhesion reinforcement triggered by the unfolding of the talin molecules that leads to the vinculin binding and in turn to an increased integrin density at an FA site [21,22].While an augmented motor-clutch theory with adhesion reinforcement successfully captures the monotonic increase of the traction force with the matrix stiffness [23], it cannot explain the non-monotonic dependence of the migration speed on the ECM rigidity measured in the experiments [5][6][7][8].This suggests that the mutual dynamics of the many FA sites across the cell along with the cytoskeletal dynamics must be considered to relate the migration speed to the traction forces and in turn to the matrix rigidity, rather than at a single FA site.To reproduce the cell migration dynamics as a coordinated sequence of peripheral protrusion and contractions [4,24,25], whole-cell level formulations with multiple FA sites have been developed, including our recent theory that integrates the chemomechanical coupling between the Rho GTPase concentrations, the FA sites, the cytoskeletal network and the nucleus dynamics [26][27][28].However, none of these studies have incorporated the adhesion reinforcement mechanism into the whole-cell dynamics.Importantly, we anticipate that a mere addition of the adhesion reinforcement effect to the whole-cell-scale model would still lead to a monotonic rigidity-speed relationship, implying further subtleties associated with the collective intracellular dynamics under adhesion reinforcement.
Here we generalize our multiscale theory in Shu & Kaplan [28] to elucidate the spreading and migration dynamics of adhesion-reinforced cells on viscoelastic substrates, overcoming the limitations of our previous model without adhesion reinforcement.Our model demonstrates how the nonlinear mechanics of the nucleus-cytoskeletal network complex must play a critical role in the experimental spreading and non-monotonic stiffness-speed profiles: large cell spreading on stiff substrates induced by the augmented local traction forces due to adhesion reinforcement, an effect accurately captured by our model, deforms the nucleus [29] and increases both the cytoplasmic viscosity and cytoskeletal stiffness [8,[30][31][32].Consequently, the drag force on the nucleus increases, making it harder for the cell to translocate the nucleus on stiff substrates, as evidenced by a theoretical non-monotonic stiffness-speed relationship in quantitative agreement with experiments.By using the standard linear solid (SLS) model for the substrate viscoelasticity, we also investigate the effect of the substrate stress relaxation on cell migration across a broad stiffness range, which is experimentally well documented [23,[33][34][35].We show that, on a substrate with high elastic stiffness, fast stress relaxation (low viscosity) may pre-empt the effect of adhesion reinforcement whereas slow relaxation (high viscosity) may promote it.That way, our model quantitatively reproduces the experimental spreading area and migration speed profiles of the HT-1080 human fibrosarcoma cells [27].Furthermore, we demonstrate that adhesion-reinforced cells persistently move up the stiffness gradients, i.e. durotaxis [36][37][38], in good agreement with the experimental migration patterns of the MDA-MB-231 cells [16].We also reveal the effect of the viscosity gradients on directed migration, a.k.a.viscotaxis, which sheds light on the migration of human mesenchymal stem cells from high to low loss moduli on collagen-coated polyacrylamide gels [39].
Our theory provides a rigorous understanding of how the nucleus drag as a function of cell spreading gives rise to the experimentally observed complex migration dynamics, in contrast with the previous phenomenological treatments that assumed a direct functional relationship between the drag force, the traction force and substrate stiffness [40,41].Furthermore, because our model takes into account additional intracellular components beyond just the FA sites, it provides a plausible mechanism for the emergent non-uniform stiffness-speed relation that complements the cell model which relates cell speed only to the resultant traction force from the collective adhesion-reinforced FA dynamics [42].

Methods
Mesenchymal migration entails the spatio-temporal coordination of front protrusion driven by actin polymerization and focal adhesion (FA) followed by rear contraction due to focal deadhesion and actomyosin activity.The front-rear symmetry is primarily broken by the chemical polarization of the active and inactive Rho GTPase proteins, which govern the dynamics of the FA sites.
Here we generalize our previous multiscale framework for the talin-low cell migration, which accounted for this complexity by considering the feedback between the biological components in figure 1a, by including the effect of focal adhesion reinforcement due to talin unfolding (figure 1b) [28].Our modified multiscale model takes into account the nonlinear strain-stiffening of the cytoskeleton and the augmented viscous drag force on the deformed cell nucleus under cell flattening, as well as uses the SLS model for the ECM viscoelasticity.At the subcellular scale, a motor-clutch model with adhesion reinforcement is employed to determine the effect of substrate rigidity on traction forces at each FA site (figure 1c).Since cells operate in low Reynolds numbers, we only consider the viscous forces arising from the nucleus-cytoplasm and cell membrane-substrate frictions at the cellular scale.In our model, the forces across the two scales are transmitted by the cytoskeleton, and the augmented viscous drag force on the deformed cell nucleus in turn determines the active motion of the cell (figure 1d,e).

Equations of motion at subcellular scale
We first present the calculation of the cell vertex displacements due to the FA site dynamics and the mechanical balance at the corresponding vertex.For a vertex i with position vector x i (t), the local spreading velocity V i s ; _ x i needs to be determined at every time step.In our model, the radial component V i s,n ; V i s Á ni is controlled by the active processes, i.e. actin polymerization and actomyosin contraction regulated by the Rho GTPase proteins, while the polar spreading speed V i s,t ; V i s Á ti is set by the balance between the polar components of the passive forces due to the cell deformations (n i , ti defined in figure 1d ).
The active radial spreading of a vertex with a speed V i s,n is fuelled mainly by F-actin formation with a polymerization speed V i p and counteracted by the retrograde G-actin flow with a speed V i r , which yields The polymerization rate at the vertex i is given by the ratio of the active Rac1 concentration R i a to its mean hR i a i averaged over all vertices as V i p ; ðR i a =hR i a iÞV 0 p , where V 0 p is the reference polymerization speed.The retrograde actin flow speed V i r is promoted by the net resistance force against protrusion F i p due to the cell membrane elasticity and myosin contractions (figure 1b,c) and impeded by an elastic restoring force F i c due to the formation of molecular bonds by proteins such as integrins, talin and vinculin between F-actin and ECM.We thus propose a phenomenological relation where V 0 is the unloaded myosin motor speed, N i m is the number of active myosin motors, and f m is the force that stalls the activity of one myosin motor.Owing to the increased iteration requirements for calculating V i r when reaching its minimum during loading-unloading focal adhesion cycles, we enforce the condition V r ¼ maxðV r , 0Þ that preserves the accuracy while reducing the computational load (see electronic supplementary material, §S1 and figure S1, S2).RhoA is known to induce myosin motor activation, leading to stress fibre formation and contractility [43,44].Therefore, we assume that the myosin motor number is controlled by the active RhoA concentration r i a and its mean hr i a i averaged over all vertices as N i m ðtÞ ; ðr i a =hr i a iÞN 0 m , where N 0 m is the reference myosin motor number.Note that this assumption guarantees the constancy of the total motor number in the cell, in accordance with the conservation of mass [16,26].In equation (2.2), the elastic restoring force F i c is determined from the FA dynamics with adhesion reinforcement, which is triggered by the stiffness of the viscoelastic substrate.On the other hand, the protrusion force F i p is set by the local force balance at each vertex in the presence of the nonlinear strain stiffening of the cytoskeleton.To calculate these two force strengths, we detail each of those processes next.

Focal adhesion dynamics with adhesion reinforcement
To account for the adhesion reinforcement due to talin unfolding, we extend the augmented motor-clutch model for FA dynamics introduced by Elosegui-Artola et al. [21] and Gong et al. [23] in order to calculate the clutch force F i c at each vertex.By denoting average displacements of all bounded clutches at the filament end by x i r ðtÞ and the displacement of the substrate by x i sub , the engaged clutch is represented by a Hookean spring with tension [45,46].The clutch force f i c displays oscillatory dynamics driven by the cyclic process of binding, loading and unbinding within the system (electronic supplementary material, figure S3).At any instant t, the unbounded clutches associate with a rate k i on .For the talin-low cells, a constant association rate k 0 on is typically assumed, corresponding to a constant clutch binding timescale t 0 on ; 1=k 0 on [17,23,28].However, a cell can activate adhesion reinforcement once it probes the elevated loading rate and high force per clutch before talin unfolding takes place.We thus introduce the adhesion reinforcement by assuming that, when the timeaveraged clutch force hf i a i t l ; Ð t l 0 fi c d t=t l (τ l : variable focal adhesion lifetime) without adhesion reinforcement is above a threshold force f cr , the clutch binding rate will increase per [21,23], where ζ is a characteristic inverse force scale.Since hf i a i t l must be obtained under the assumption that the adhesion reinforcement be absent, we separately calculate it by solving an isolated motor-clutch system with ζ = 0 (equations (2.2)-(2.6),see inset in electronic supplementary material, figure S4 and  §S2), considering the updated protrusion force F i p and the updated number of clutches N i c and motors N i m at vertex i.The clutch force fi c and the time t are variables in isolated motor-clutch modelling, distinct from the variables in our whole cell modelling.Owing to the time dependence of F i p , N i c and N i m , the variable hf i a i t l at every vertex must be updated at every N time steps.Our simulations have demonstrated that using N = 4000 time steps provides accurate predictions on cell migration speeds on stiff substrates, comparable to results obtained from more frequent updates of hf i a i t l (see electronic supplementary material, figure S5), while keeping computational costs reasonable.
Once formed, the molecular clutches must unbind at a dissociation rate k i off that depends on the clutch tension f i c .To that ; D is the diffusion constant of the membranebound species (equation (2.15)).The angle between the two membrane sections at a vertex is denoted by θ i .(e) The principal dimensions l max , l med , l min of a deformed nucleus.
end, we use the functional form Here, k r0 and k c0 denote the unloaded off-rate and the unloaded catch-rate, respectively, f r0 is the characteristic rupture force, and f c0 is the characteristic catch force.It follows that the fraction of the engaged clutches (0 ≤ P i (t) ≤ 1) is governed by the mean-field rate equation [23,45,46], Denoting the number of available clutches as N i c , the total clutch force at the vertex i is then given by F i c ¼ P i N i c f i c .Here we incorporate the critical role of the Rac1 proteins in focal complex assembly by relating N i c to the local Rac1 concentration, i.e.N i c ðtÞ ; ðR i a =hR i a iÞN 0 c , where N 0 c denotes the reference clutch number [47][48][49].This assumption implies that the total count of clutches P N i¼1 N i c remains constant at any given moment, thereby adhering to the conservation of mass [16,26].The mechanical equilibrium condition at the cell-substrate interface demands that the total force sustained by the engaged clutches must be balanced by the substrate deformation, leading to In equation (2.5), the substrate displacement x i sub is an unknown to be determined from a constitutive model for the substrate viscoelasticity, which we focus on next.

Constitutive model for substrate viscoelasticity
Our previous implementation used the classical Kelvin-Voigt model for the substrate relaxation dynamics [28].However, since the Kelvin-Voigt model predicts a very rigid non-physical behaviour when t < τ r (τ r : substrate relaxation timescale) [50], it can result in a premature adhesion reinforcement.Comparatively, the SLS model allows for a better physical modelling by including a Maxwell arm (figure 1) [51].It has also been demonstrated that the SLS model can capture the prominent relaxation timescale of the viscoelastic substrates fabricated by combining covalent and supramolecular cross-linking [23].The SLS model expresses the constitutive relationship between x i sub and the substrate force F i sub as (figure 1c) where K e is the elastic stiffness at t → ∞, γ is the substrate viscosity, and K a is the additional stiffness that governs the substrate relaxation with a timescale τ r ≡ γ/K a .A stress relaxation test from a constant strain yields the instantaneous and long-term stiffness of the substrate as K t→0 = K e + K a and K t→∞ = K e , respectively.The instantaneous stiffness K t→0 characterizes the initial elastic response of the substrate when the viscous deformation and stress relaxation have not yet taken place in the limit t → 0. Thus, cells with a very short focal adhesion lifetime (τ l < τ r ) can only sense the instantaneous stiffness.By contrast, the longterm stiffness K t→∞ refers to the residual substrate stiffness after the viscous stresses have relaxed.

Cytoskeletal stiffening
The passively deforming cell cytoskeleton, which consists of microtubules and intermediate filaments, is represented by multiple springs in our model (figure 1d).Here we will assume that these cytoskeletal 'springs' undergo strain-stiffening during large spreading events and thus exhibit nonlinear elasticity.This assumption is backed by the experiments performed, e.g. on NIH-3T3 fibroblasts, that reveal a strong correlation between the cell rigidity and cell area during large spreading events [8].Assuming this behaviour is mechanically driven and thus must be generic across animal cells, we propose a phenomenological equation for the cytoskeletal stiffness K cs in terms of the cell area A: denoting the position vector of the cell nucleus by x nuc and defining the length of a cytoskeletal spring as r i ≡ |(x nuc − x i )|, the cytoskeletal restoring force in our model is governed by the differential equation where K b cs denotes the baseline stiffness.The term ΔK cs e βA is introduced to describe the exponential hyperelastic behaviour [52,53], with the values of ΔK cs and β obtained through linear regression of the experimental data in [8] (electronic supplementary material, §S3 and figure S6A).

Local mechanical equilibrium
We enforce local force balance at each vertex to determine the radial protrusion force F i p and the polar spreading speed V i s,t .The cell membrane possesses a notable mechanical rigidity that enables it to endure a variety of stresses, which is vital for maintaining the integrity of cell structures [41,[54][55][56].We model it as a closed loop of Hookean springs with stiffness K m between neighbouring vertices.The tensile forces acting on the ith vertex can be expressed as , where l i þ and l i À are the distance vectors between vertex i and its two neighbours, and l 0 is the initial distance between two adjacent vertices.The drag force on vertex i can be expressed as defines the average membrane length about the ith vertex.The drag force can originate, e.g. from the hydraulic resistance of the extracellular medium and the viscous response of the cortex beneath the cell membrane.Thus, the value of the drag coefficient h i m in our model is within the range of the approximate hydraulic resistance coefficient and the range of the viscous coefficient used in other twodimensional models (electronic supplementary material, §S4, figure S7) [54,56,57].In addition to the membrane forces, each vertex experiences a protrusion force F i p and a cytoskeletal radial force F i cs .Thus, the net force balance at each vertex in the radial direction ni and the polar direction ti is given by and Equations (2.8a) and (2.8b) yield F i p and V i s,t , respectively.Altogether, equations (2.1)-(2.8b)fully determine the vertex spreading velocities V i s when the nucleus displacement x nuc and the dynamical GTPase concentrations at each vertex are computed at the cellular scale.

Equations of motion at cellular scale
Next, we explain the global mechanical equilibrium that governs the nucleus motion and the intracellular Rho-GTPase dynamics, which directly influence the aforementioned subcellular processes.

Viscous drag on deformed nucleus
Our model quantifies cell translocation by the net translation of the cell nucleus, which balances the forces between the cytoskeletal microtubules and intermediate filament bundles.At the cell periphery, these cytoskeletal complexes are linked to the F-actin at the FA sites, transmitting the net traction force from the ECM (F i sub ) and the protrusion force from the cell boundary (F i p ) to the cell nucleus (figure 1).This tight linkage to the rest of the cell can deform the nucleus when the cell flattens under large spreading on stiff substrates [29].To balance these peripheral and cytoskeletal forces, the nucleus undergoes viscous drag within the cytoplasm while it deforms under cell flattening.The viscous drag force on a particle is commonly described by the product of a drag coefficient, a characteristic particle size and the velocity of the particle relative to the surrounding medium.Note that the cell peripheral dynamics cannot induce a net flow into the cell interior since the viscoelastic relaxation timescale of the cytoplasm τ cp ∼ 0.1−1 s is much smaller than the cell migration timescale T ≡ L/V 0 ≈ 260s (L = 2πr 0 : size of a circular cell, V 0 : retrograde flow speed; electronic supplementary material, table S1) [58].As a result, flow around the nucleus should predominantly be induced by the dynamics of the nucleus and other organelles, as well as the cytoskeleton.In the absence of a detailed characterization of the organelles and the cytoskeleton, we assume that the drag force is only proportional to the nucleus speed _ x nuc [28,41,54] where η cp represents the viscosity of the cytoplasm and L nuc is the characteristic particle size.For simplicity, we take L nuc ; 6pr 0 nuc corresponding to Stokes' flow as a first-order approximation to the cytoplasmic domain within finite cell height (r 0 nuc : the initial radius of the nucleus).The shape factor f shape quantifies the deviation of the nucleus from a perfect sphere with f shape = 1.The Corey shape function can be used to determine f shape based on the nucleus's three principal lengths (figure 1e) [59,60] where the exponent α = 0.09 was obtained by fitting the experimental drag coefficient of non-spherical particles under Stokes flow [60].For the deformed nucleus, the aspect ratio defined by the longest and the shortest dimensions can be related to the cell spreading area by δ ≡ l max /l min = χA + 1, where the exponent χ = 0.0024 is obtained by fitting the experimental cell shape data in Li et al. [61] (electronic supplementary material, figure S6B).The intermediate dimension l med corresponds to the length of the minor axis of the nucleus on the x-y plane, and it was found to have a constant ratio to the length of the major axis in experiments (l med = 0.8l max ) [29].Consequently, the nucleus shape factor can be simply related to the cell spreading area A by Experiments have demonstrated that the cytoplasm viscosity, like the strain stiffening, is strongly correlated with the cell spreading area [31,32].These experiments also indicate that the cell viscosity and cell stiffness exhibit the same trend with increasing substrate stiffness.Therefore, we assume an area-dependent cytoplasm viscosity (similar to that in equation (2.7)) as

:12Þ
where h 0 cp is the cytoplasm viscosity on soft substrates [62] and Δη cp controls the rate of viscosity increase with the cell spreading area A. The fitting parameter Δη cp ensures that the maximum cytoplasmic viscosity is approximately twice its value on soft substrates as indicated by experiments [31].The values of the Δη cp may vary depending on the cell type and physiological conditions, which could explain the quantitative differences in non-monotonic stiffness-speed relationships among different cell types [5,6,29].
At the frame of the nucleus, the condition for the mechanical equilibrium between equation (2.9), the cytoskeletal and the subcellular forces yields the nucleus migration velocity, _ x nuc , equivalent to the cell velocity in our model.

Reaction-diffusion dynamics of GTPase concentrations
As with our previous work [28], here we adopt the biochemical reaction-diffusion equations introduced in [63,64] to describe the dynamics of active and inactive GTPases.Since the active forms of the GTPases are predominantly associated with the cell membrane, which also serves as a major site for the conversion between active and inactive forms, we track the volume fractions of the signalling proteins in three forms [65,66]: the active membrane-bound form ðG i a ðtÞ ; fR i a ðtÞ, r i a ðtÞgÞ, the inactive membrane-bound form ðG i in ðtÞ ; fR i in ðtÞ, r i in ðtÞgÞ, and the inactive cytosolic form (G cp (t) ≡ {R cp (t), ρ cp (t)}).Given the rapid diffusion of inactive proteins in the cytosol, we assume G cp (t) remains uniformly distributed in the cytosol at all times.Both active and inactive membrane-bound forms diffuse with a diffusion constant D across the vertices.The corresponding diffusive fluxes are given by Fick's Law in a finite difference formulation as Equation (2.14) takes into account the effect of the deformed cell shape on the diffusive flux by updating the vertex coordinates at each time step.Different forms of the proteins on a vertex are interconvertible with the inactive-to-active rates A i G , active-to-inactive rates I i G , and inactive-to-cytosolic association and disassociation rates M þ G , M À G (figure 1d).Altogether, the reaction-diffusion kinetics is governed by ð2:15Þ The total amounts of Rac1 and RhoA are constant due to the conservation law The active Rac1 and RhoA GTPase volume fractions R i a , r i a regulate the cell migration by controlling the actin polymerization speed V i p (equation (2.1)), the number of myosin motors N i m (equation (2.2)) and the clutch number N i c at each FA site (equation (2.5)).Different from Merchant et al. [63] and Zmurchok & Holmes [64], the rate terms A i G in our formulation account for the reverse coupling from the cell deformation to the signalling pathways in addition to the mutual inhibition of the Rac1 and RhoA (appendix A.1). Since the filopodial protrusions of a mesenchymal cell grow and shrink at timescales comparable to the migration times, this chemo-mechanical feedback prevents the Rac1 and RhoA dynamics from reaching a steady bistable polarized state.Instead, when the transient chemical polarity ceases, our simulation algorithm reinstates the polarization stochastically to sustain the random migration patterns in many mesenchymal phenotypes [26,27].

Simulation procedure
We represent the initial vertex configuration of the cell as a regular 16-sided polygon based on the mesh convergence study in [28], which suggests that using N = 16-24 ensures decent numerical accuracy.All initial forces, velocities, displacements and substrate deformations are set to zero.Since chemical signalling drives cell polarization followed by directional locomotion on uniform substrates [28], a non-uniform initial distribution of the active Rac1 protein volume fraction R i a is enforced with a higher value at the cell front.Likewise, a polarized initial distribution of the active RhoA protein volume fractions r i a is taken royalsocietypublishing.org/journal/rsifJ. R. Soc.Interface 20: 20230317 with accumulation at the cell rear (manual polar symmetry breaking).In contrast, the initial conditions for R i in , r i in are uniform at time t = 0, and the cytoplasmic form R cp , ρ cp can be determined from the total volume fraction conservation (electronic supplementary material, figure S9).For directed migration on non-uniform substrates, uniform initial signalling distributions must be assumed since the polarity will be set by the gradients of the ECM stiffness or viscosity.Starting from these initial conditions, the vertex spreading and the nucleus velocities are calculated by solving equations (2.1)-(2.13),and the Rho-GTPase volume fractions are updated by solving equations (2.14), (2.15), (A 1) and (A 2) by following the algorithm in electronic supplementary material, figure S4.When the cytoplasmic inactive Rac1 volume fraction R cp reaches a steady state, we reinitialize R i a and r i a stochastically to mimic the random nature of the protrusion formations in the mesenchymal cells.To incorporate the adhesion reinforcement (equation (2.3)), we fix ζ = 0.5 pN −1 since it leads to a good agreement between the simulated spreading areas and the experimental data in Solon et al. [8] (electronic supplementary material, figure S10).The parameters ΔK cs = 2.8 pN μm −1 and β = 0.00185 (equation (2.7)) and χ = 2.4 × 10 −3 (equation (2.11)) provide the best fit to the experimental data from Solon et al. [8] and are kept constant in all simulations (electronic supplementary material, figure S6).All fixed simulation parameters are listed in electronic supplementary material, table S1.
We present the simulation results in dimensionless form by introducing a position scale L ≡ 2πr 0 = 10π μm, retrograde flow speed scale V 0 = 120 nm s −1 , migration timescale T ≡ L/V 0 ≈ 260 s and force scale F 0 ; N 0 m f m ¼ 200 pN.For each simulation, we run M > 10 6 time steps with a time step of approximately Δt = 1.5 × 10 −3 s (Δt ≈ 5.7 × 10 −6 in dimensionless units), corresponding to the migration dynamics over at least an hour in real units.The spreading area A is calculated as the average cell area over the last 10 6 time steps in a simulation.We define the migration speed V as the nucleus trajectory length X ≡ |x nuc (MΔt) − x nuc (0)| divided by the total time t total ≡ MΔt, i.e.V ≡ X/t total .The ranges of the ECM material parameters used in the simulations are listed in table 1.For each set of the ECM parameters, we run n simulations to compute the sample mean spreading area, A, and the mean migration speed, V.The effects of the ECM viscoelasticity on the spreading area, migration speed and migration direction are validated by the experimental data taken from [6,16,27].

Simulations reproduce experimental cell spreading and migration speed dependence on substrate stiffness
With increasing stiffness K e on an elastic substrate (γ = 0, K a = 0), the cell spreading area increases monotonically, in agreement with the experiments on the adhesion-reinforced U373-MG human glioma cell spreading on polyacrylamide hydrogels (figure 2a) [6].This can be understood by inspecting the FA dynamics at the subcellular scale: according to equations (2.2), (2.5), the cell spreading speed at a vertex V i s,n is set by the competition between the polymerization speed V i p and the actin retrograde flow speed V i r , which is counteracted by the stiffness-dependent focal traction force (appendix A.2) ) is derived for an isolated motor-clutch system where the engaged clutch fraction is assumed to be saturated on a compliant substrate (P i ≈ 1; K e ( N 0 c K c ).For sufficiently low K e where the adhesion reinforcement mechanism is not triggered (see equation (2.3)), the traction force starts building up at an FA site after the molecular clutches form at a timescale t 0 on ; 1=k 0 on , and develops with a characteristic time τ l (equation (3.1)).The comparison between the two timescales ðt 0 on ¼ t l Þ defines a threshold stiffness K 0 ; N 0 m f m k 0 on =V 0 % 5 pN nm À1 (table 2).On soft substrates (K e ≤ K 0 ), since the force build-up time τ l is much longer than the clutch formation time τ on , the adhesion force mostly remains low, i.e.F i sub , N 0 m f m (equation (3.1); figure 2b).For a stiff substrate with K e > K 0 , the decrease in the force build-up time τ l must cause a steep increase in the average force per bond (hf a i t l .f cr ), triggering the adhesion reinforcement mechanism with an augmented binding rate k on (equation (2.3), electronic supplementary material, figure S3).The reinforced binding rate can even saturate the bounded clutches on a very stiff substrate, i.e.P i → 1 (electronic supplementary material, figure S3F).Moreover, a higher hf a i t l on stiff substrates leads to higher disassociation rates k off , which must balance the binding rate in equation (2.4) to reach a steady bound clutch fraction P i ∼ 1.This condition yields an estimate for the instantaneous force per bond as f i c / zðhf a i t l À f cr Þf r0 (appendix A.3).This allows us to estimate an upper bound of the substrate traction force on a stiff substrate as in agreement with the simulations (figure 2b).This elevated traction at K e > K 0 across the cell is the main driving factor behind enhanced cell spreading.Furthermore, the Rho GTPase polarization leads to a higher traction force at the cell front than at the rear since the motor-clutch number N i c is controlled by the Rac1 concentration (figure 2b).Owing to the chemically induced polarization, the cell acquires a larger net traction force on stiff substrates than on soft ones (figure 2c).
The cell motility is governed by the competition between the net traction force and the viscous drag of the nucleus since the remaining forces are negligible in equation (2.13) (electronic supplementary material, figure S13).Because large cell spreading on stiff substrates can deform the nucleus (figure 2d) and alter the cytoplasmic viscosity (equation (2.12)), the viscous drag force on the nucleus must increase, resisting the net traction force.Consequently, the migration efficiency decreases beyond the threshold stiffness K 0 , creating a non-monotonic stiffness-speed profile, which also quantitatively reproduces the U373-MG human glioma cell speed on polyacrylamide hydrogels (figure 2e) [6].In our simulations, ignoring the nucleus shape change and assuming constant cytoplasmic viscosity recovers a monotonic relationship between the migration speed and substrate stiffness, which validates their role in the nonmonotonic speed-stiffness relation (figure 2f ).Furthermore, we have performed parametric studies to explore the relative impact of the nucleus shape change and the cytoplasmic viscosity increase on the non-monotonic speed-rigidity relationship (electronic supplementary material, figure S14).We found that while both factors influence the migration speed on soft substrates, the effect of cytoplasmic viscosity increase becomes more significant on stiff substrates, emphasizing its critical role in capturing the non-monotonic relationship.We also examined the potential impact of cytoplasmic flows on nucleus motion by modifying equation (2.9) to consider the centroid velocity.Electronic supplementary material, figure S8 demonstrates that this approximation does not substantially influence the efficiency of cell movement over long times.Overall, because glioma cells can migrate with or without adhesion reinforcement [16,69], our current results complement our previous work in [28] by explaining how a non-monotonic speed-stiffness relation emerges while the traction forces remain monotonic.

Substrate viscosity affects cell spreading area and migration speed
We next discuss the effect of substrate viscosity on the cell area and migration speed on a substrate that is soft in the long time limit, K e = 0.1 pN nm −1 .When K e + K a < K 0 (soft regime), the adhesion reinforcement effect can be neglected even on viscous substrates.However, when K e + K a ≥ K 0 (stiff regime), a large viscosity may trigger the adhesion The average dimensionless substrate traction force versus K e /K 0 at vertex 0 (cell front) and vertex 8 (cell rear).The red dash-dotted line gives the upper bound of the substrate traction force on very stiff substrates.The calculation of the averaged substrate traction forces hF i sub i t total from the full focal adhesion cycles is illustrated in electronic supplementary material, figure S11.(c) The dimensionless mean and standard deviations of the net traction force versus K e /K 0 .The calculation of the time-averaged net traction force is illustrated in electronic supplementary material, figure S12.(d ) Cell-spreading induced nucleus viscous drag increase f shape h cp =h 0 cp versus K e /K 0 , which is calculated by using equations (2.9), (2.11), (2.12) for the mean spreading area values in a. (e) The simulated dimensionless mean migration speed V=V 0 versus K e /K 0 (full curve) and the experimental data for the U373-MG human glioma cells (dots and error bars) [6].( f ) The influence of the shape factor f shape of the deformed nucleus and the area-dependent cytoplasm viscosity η cp on the migration speed (equations (2.11), (2.12)).The mean values or the standard deviations are calculated over n = 10 simulations at each data point in a, c-f.Table 2. Characteristic scales for spreading and migration dynamics.For a viscoelastic substrate, K a is the additional stiffness, and γ is the viscosity in the SLS model.Additionally, N 0 m is the reference myosin motor number, V 0 is the unloaded myosin motor speed, f m is the force that stalls the activity of one myosin motor, and k 0 on denotes the reference clutch binding rate.
quantity clutch binding timescale a clutch lifetime a threshold stiffness substrate relaxation timescale optimal viscosity definition t 0 on ; 1=k 0 on a The scales are valid for substrates with a stiffness lower than K 0 .
royalsocietypublishing.org/journal/rsif J. R. Soc.Interface 20: 20230317 reinforcement regime and alter the spreading behaviour.We discuss both regimes below.On substrates with a low instantaneous stiffness (K e + K a < K 0 ), the cell spreading area varies in a non-monotonic manner with increasing viscosity (figure 3a,b).This can be explained by examining the FA dynamics at the subcellular scale, similar to the cell spreading on elastic substrates.At an FA site without adhesion reinforcement, time averaging the traction force (see electronic supplementary material, figure S15A,B) over a clutch lifetime τ l leads to (appendix A.2) where the timescales t l , t 0 on , t r are defined in table 2. Equation (3.2a) indicates that the average traction force at an FA site must increase with increasing viscosity when t r t 0 on .At t r .t 0 on , since slower substrate relaxation leads to the FA lifetime reduction when (electronic supplementary material, figure S15C) the average traction force should exhibit a decreasing profile with a lower bound that is estimated in equation (3.2b).That way, equation (3.2) indicates the presence of a non-monotonic traction force-viscosity profile in an isolated motor-clutch system on soft substrates.Note that when t 0 on , t r , t l , the average substrate force is not amenable to an analytical approximation.Still, given the biphasic profile suggested by equation (3.2), it is plausible to assume that the average substrate force reaches a maximum when t r % t 0 on , equivalently at an optimal viscosity g 0 ; K a t 0 on (electronic supplementary material, figure S15C), thereby maximizing the cell spreading area.To further validate our model, we simulated spreading areas of human fibrosarcoma cells HT-1080 on viscoelastic substrates consisting of interpenetrating alginate networks (IPN) [27].The experimental material parameters of the substrates with different relaxation rates were obtained by data fitting to the stress relaxation tests in [27] (electronic supplementary material, §S3 and figure S16).Clearly, our whole-cell model accurately replicates the experimentally observed spreading areas of the HT-1080 cells on soft substrates with different stress relaxation rates (figure 3c), indicating that cell spreading is suppressed with increasing viscosity when t r .t 0 on .At a larger instantaneous substrate stiffness K e + K a ≥ K 0 , the cells only perceive the long-term stiffness K e and maintain small spreading areas at a low substrate viscosity γ (figure 3d).For higher γ, the slowly relaxing substrate maintains the stiffness K e + K a ≥ K 0 at later times, sustaining the adhesion reinforcement regime and thus promoting cell spreading.Moreover, because the adhesion reinforcement effect becomes more pronounced on stiffer substrates, a larger additional stiffness, such as K a = 20 pN nm −1 further increases the cell spreading area at high γ (figure 3d).
Figure 3.Effect of substrate viscosity on spreading area and migration speed.(a) Contour plot of the mean dimensionless cell area A=L 2 on viscoelastic substrates as a function of the additional stiffness K a and substrate viscosity γ (K e = 0.1 pN nm −1 ).The white line denotes the threshold stiffness K 0 .(b) A=L 2 and its standard deviation (shaded area) versus the ratio of the material relaxation timescale τ r ≡ γ/K a to the clutch binding timescale t 0 on ; 1=k 0 on on soft substrates.(c) For a fixed instantaneous stiffness K e + K a = 2.6 pN nm −1 , the simulated spreading area A and its standard deviation (n > 50) and experimental data for HT-1080 cells (data from [27]).By fitting the stress-time data in the stress relaxation test, the viscosities of the fast, intermediate (med) and slow relaxing substrates are found as γ fast = 1 pN s nm −1 , γ med = 60 pN s nm −1 , γ slow = 360 pN s nm −1 (electronic supplementary material, figure S16).(d ) The dimensionless spreading area A=L 2 and its standard deviation as a function of γ for different K a > K 0 .Since the clutch binding timescale τ on depends on the adhesion reinforcement on stiff substrates, A=L 2 is plotted against the substrate viscosity γ.(e) The mean migration speed V and its standard deviation in the simulations (n > 50) versus experiments of HT-1080 cells (data from [27]) on fast-relaxing and slow-relaxing soft substrates, for which the material parameters are extracted from fitting the stress-time data (electronic supplementary material, figure S16).( f ) The dimensionless mean migration speed V=V 0 and its standard deviation as a function of γ for two additional stiffness values K a = 8 pN nm −1 and K a = 20 pN nm −1 .The contour plot in (a) contains 180 data points.The mean and standard deviations are calculated with n ≥ 5 simulations at each point in b, d, f.In c, e, statistical analyses were conducted using the t-test to compare simulation results with experimental data.Levels of statistical significance ( p) are indicated as follows: # p , 0:0001 and ***p < 0.001.
royalsocietypublishing.org/journal/rsif J. R. Soc.Interface 20: 20230317 As opposed to the spreading area that can be understood at the single FA level, the change in the cell migration speed as a function of γ must be examined by virtue of the net traction force P i F i sub at the cellular scale.At an instantaneous substrate stiffness K a + K e < K 0 , where the adhesion reinforcement is not activated, the high viscosity of γ > γ 0 can result in a decrease of the FA lifetime and thus a lower traction force (appendix A.2 and electronic supplementary material, figure S15D).Our simulations reveal that migration on fast relaxing substrates (t r t 0 on ) is more efficient than on slow relaxing substrates (t r ) t 0 on ) (figure 3e).This is in quantitative agreement with the migration speeds of the MDA-MB-231 human breast cancer adenocarcinoma cells on fast and slow relaxing IPN [27].A more comprehensive parametric study demonstrates the presence of a biphasic speed-viscosity relation (electronic supplementary material, figure S17A,B).The optimal viscosity corresponding to the maximum migration speed tends to be slightly higher than γ 0 (for maximum cell spreading), primarily due to the prolonged FA lifetime at the front of the cell (appendix A.2). On substrates with high instantaneous stiffness (K a + K e ≥ K 0 ), the net traction force increases with viscosity due to the adhesion reinforcement effect (electronic supplementary material, figure S17C).Since cells have limited spreading areas due to the weak engagement of the adhesion reinforcement regime at K a = 8 pN nm −1 (figure 3d), the change in the nucleus drag force is marginal.Thus, we obtain a nearly monotonic increase in the migration speed with viscosity (figure 3f ).However, with a very large additional stiffness of K a ≫ K 0 , strong adhesion reinforcement on the viscous substrate results in very large spreading areas (figure 3d), which significantly increases the nucleus drag forces and in turn leads to a biphasic speed-viscosity profile (figure 3f ).

Migration direction is set by stiffness or viscosity variation on non-uniform substrates
Having elucidated the effect of substrate stiffness and viscosity on cell area and migration speed, we investigate the sensitivity of cell migration to a stiffness gradient on elastic substrates or a viscosity gradient on viscoelastic substrates.To eliminate the effect of chemical polarity on migration direction in the simulations, we specify the initial conditions of a chemical apolar cell by assigning a random distribution of the active Rac1 concentration R i a and setting the other membrane-bound Rho GTPase concentrations r i a , R i in , r i in constant (electronic supplementary material, figure S18).Migrating cells exhibit limited migration speeds (V , 50 mm h À1 ), leading to a substantial amount of time required for simulating cell locomotion over long distances.Therefore, for the sake of computational efficiency, we adopt a 200 × 300 μm domain at about a 1 : 10 scale of the experimental platform (figure 4a).We compared our migration simulations of 83 individual cells on elastic gradient substrates with the experimental migration patterns of the MDA-MB-231 cells on polyacrylamide hydrogels with a stiffness gradient 0.5-22.0kPa [16].We record the migration trajectory of each simulated cell for 7.2 h, which is 1/10 of the experimental observation time (72 h).Owing to the difficulty in quantifying the substrate gradient in the experiments, we assume a constant gradient rK e % 72 pN nm −2 in the +y-direction in the simulations; non-dimensionalization by the initial cell diameter 2r 0 and the threshold stiffness K 0 yields the unitless gradient r Ke % 0:14 , 1.The simulated cell trajectories indicate that all cells with random locations at t = 0 move up the stiffness gradient regardless of their initial positions, i.e. the cells undergo robust durotaxis (figure 4a).Also, the percentages of the simulated cells in three substrate portions demonstrate the strong migration trend towards the stiffest region (figure 4b), in good agreement with the gradient sensitivity of the MDA-MB-231 cells with talin-and vinculinmediated FA formations [16].We further confirm that almost all cells eventually migrate to the stiffest region after more than 12 h (figure 4b).The robust durotaxis behaviour is mainly due to the monotonic increase of the traction force as a function of the substrate stiffness at an FA site due to the adhesion reinforcement (figure 2b).In other words, the FA sites attached to the stiffer regions generate a higher traction force than those attached to the softer regions on a non-uniform substrate, driving the cell up the stiffness gradient.
We next sought to understand the impact of substrate stress relaxation on the directed migration of cells on soft substrates.In our simulations, the two stiffness values in the SLS model are set to be K a = 2.9 pN nm −1 and K e = 0.1 pN nm −1 .We assume a non-uniform viscosity γ < γ 0 ≈ 1 pN s nm −1 with a constant gradient rg ¼ 12:5 pN s m À2 in +x-direction, equivalent to rg % 0:125 in the unitless form (non-dimensionalized by γ 0 and 2r 0 ).For computational efficiency, we used a 80 × 40 μm domain, which still much larger than a migrating cell size in ±x-direction.The simulations demonstrate that cells migrate from the soft elastic region (γ → 0) to the viscoelastic region, which is also evidenced by the biased angular displacement towards the +x-direction (figure 5a,b).The same trend was observed in human mesenchymal stem cells that migrate towards regions with a larger relaxation modulus on a collagen-coated polyacrylamide gel.This behaviour is referred to as viscotaxis [39].
Once the substrate viscosity surpasses the threshold (γ > γ 0 ), cells migrate against the viscosity gradient with the angular displacements biased toward faster relaxation regions, which can be referred to as 'anti-viscotaxis' (figure 5c,d).As a control, the random cell trajectories and unbiased angular displacements were observed on a uniform viscoelastic substrate (figure 5e,f ).These findings suggest that a viscosity gradient draws the cells to a region with t r % t 0 on .The underlying mechanism for the dependence of the migration direction on the viscosity gradient can be explained by inspecting the biphasic viscosity-traction force relationship at an FA site.On soft substrates with a relaxation timescale t r , t 0 on , an increasing viscosity contributes to an enhanced substrate traction force that counteracts retrograde actin flow.This leads to a substrate-induced polarity along the gradient.By contrast, slower relaxing substrates (t r .t 0 on ) lead to decreasing clutch lifetime (electronic supplementary material, figure S15C) and the reduced mean traction forces (equation (3.2)).As a result, the edges of the cell attached to the faster relaxing substrate regions gain more traction, navigating the cell against the viscosity gradient (figure 5g).

Discussion
Our multiscale theory involves a detailed description of the adhesion reinforcement regime at each FA site and includes the cytoskeleton strain-stiffening effect, which quantitatively reproduces the large cell spreading areas on stiff substrates.Importantly, by incorporating the viscous drag force increase due to the nucleus shape change and nonlinear cytoplasmic viscosity as a function of the spreading area, we have reproduced the non-monotonic dependence of the migration speed on substrate stiffness.Our findings indicate that the interplay among larger spreading areas on stiff substrates, variable cytoplasmic viscosity, and nucleus flattening can influence the drag forces in cell migration and give rise to the nonmonotonic speed profile.Therefore, the mechanism underlying the non-monotonic speed-stiffness relation in this study is distinct from the biphasic speed-stiffness relation discussed in our previous work [28].The biphasic speed profile in [28] is attributed to the biphasic dependence of the focal adhesion force on substrate stiffness in cells that lack an adhesion reinforcement.A limitation of our previous model was that it cannot reproduce the monotonic increase of spreading area with stiffness.By contrast, our current model effectively reconciles both the observed monotonic spreading area-stiffness relationship and the non-monotonic speed-stiffness relationship, as reported in many experiments [5][6][7].Another important addition in our current study compared with our previous model is the inclusion of the cytoskeletal strain-stiffening effect.This effect plays a crucial role in accurately quantifying both cell spreading and migration efficiency of cells with adhesion reinforcement (see electronic supplementary material, figure S19).Ignoring this effect could lead to a significant overestimation of the spreading area on stiff substrates, leading to a migration behaviour that nearly comes to a halt.
Our simulations on viscoelastic substrates also reveal the influence of substrate viscosity on cell migration by tuning the engagement of the adhesion reinforcement regime.On a substrate with low instantaneous stiffness (K a + K e < K 0 ) where the adhesion reinforcement regime is absent, we have analytically shown the presence of a biphasic traction forceviscosity relationship at an FA site, with the traction force maximized on fast-relaxing substrates (t r % t 0 on ).The analytical relation helps us explain the simulation results, which quantitatively agree with the experiments.On substrates with a large instantaneous stiffness (K a + K e ≥ K 0 ), the substrate viscosity can greatly enhance the cell spreading area due to the engagement of the adhesion reinforcement regime.Notably, although royalsocietypublishing.org/journal/rsif J. R. Soc.Interface 20: 20230317 the net traction force also increases with viscosity, we observe a biphasic speed-viscosity profile with a very large additional stiffness of K a ≫ K 0 .This is because high viscosity on a very stiff substrate triggers strong adhesion reinforcement, which induces a large cell spreading area, in turn raising the nucleus drag force that counteracts the net traction force.Finally, we have investigated the influence of substrate viscoelasticity on directed migration.Our model successfully captures the robust durotaxis behaviour of the MDA-MB-231 cells with talin-/vinculin-mediated clutch reinforcement on elastic gradient substrates.These results complement our previous findings regarding the talin-low cells that operate without adhesion reinforcement and thus exhibit the cooccurrence of durotaxis and anti-durotaxis on gradient substrates.As a result, our investigations have successfully replicated the complete spectrum of directed migration behaviours displayed by MDA-MB-231 cells reported in [16].This underscores the pivotal influence of the adhesion reinforcement regime in governing cell durotactic responses on gradient substrates.Furthermore, our simulated cells also exhibit both viscotaxis and anti-viscotaxis, with all cells migrating towards the fastest-relaxing substrate region where the traction force at an FA site is maximum when t r % t 0 on .This reveals the critical role of stress relaxation in modulating the migration direction.
Our framework has inherent physical and mathematical limitations in that it is built by coupling cellular processes that are not yet fully understood.For instance, the model only takes into account the active regimes of the actin cytoskeleton at the sub-cellular level using the motor-clutch model.A rigorous modelling of the actin cytoskeleton could involve adopting the active gel theory, which characterizes dynamic structures formed by the interactions of actin filaments and associated proteins [70][71][72].Moreover, the interplay between the cytoplasmic flow and the actin cytoskeleton plays a pivotal role in cellular function and behaviour [73].While we show the minimal effect of cytoplasmic flows on the cell migration speed by approximating the net cytoplasm velocity using the centroid velocity (electronic supplementary material, figure S8), it is important to clarify that this approximation does not physically represent the cytoplasmic flow that is governed by the dynamics of cytoskeletal elements and the resulting fluid-structure interactions.Incorporating fluid mechanics into the active gel theory of the cytoskeleton can also enrich our comprehension of the dynamic nature of cells and enable us to understand certain unresolved phenomena, such as the transition between amoeboid and mesenchymal migration [74,75].
These limitations notwithstanding, our theory provides a valuable tool to understand how the intricate interplay between cellular signalling and mechanics along with cell-ECM interactions regulates migration on viscoelastic substrates.Hence, the theory paves the way for the development of novel experimental platforms for the manipulation of cell migration patterns.
Data accessibility.The Python source codes of the simulations can be accessed from the Dryad digital repository: https://doi.org/10.5061/dryad.h9w0vt4qr.
Supplementary material is available online [76].
Declaration of AI use.We have not used AI-assisted technologies in creating this article.
where K þ b governs the coupling strength between the local contractile dynamics and active RhoA concentration at vertex i.We take the value of the K þ b to be equal to the coefficients α R and β R .We assume that the active rate of Rac1 soars once the increase of the vertex angle exceeds the threshold c 0 θ 0 due to the cell contraction.Coupling Rac1 and RhoA to the mechanical deformations enables recurrent polarization and the recovery of cell deformations.Specifically, an excessive increase in θ i at the retracting cell rear builds up R i a to increase actin polymerization and thus reverses the retraction at the end of each migration step.A large decrease in θ i at the protruding cell front boosts ρ a that increases actomyosin contraction to resist further protrusion.This two-way chemomechanical coupling in our model is essential for the successive cell polarization and apolarization events while exerting minimal influence on average cell spreading area and migration speed.

A.2. Analytical derivation of the soft substrate force and deformation
We have derived an analytical expression for the substrate force F sub at a single FA site to explain the influence of viscosity on a soft substrate.The derivation focuses on isolated motor-clutch dynamics that can also be numerically solved by equations (2.2)-(2.6)(electronic supplementary material, figure S15A).For a steady-state solution, we neglect the time-dependence of the variables N c , N m , F p by setting N c ¼ N 0 c , N m ¼ N 0 m , F p = 0, which are otherwise determined by the chemo-mechanical dynamics at the cellular scale.It is important to note that our goal is not to replicate the substrate force at each FA site of the proposed whole-cell model but to decipher the mechanism of the non-monotonic profile between the viscosity and substrate force on soft substrates.Such motor-clutch analyses that ignore the details of the global dynamics prove to be effective in understanding the influence of matrix mechanics on traction forces [17,23,45,46].To that end, the analytical derivations presented are aimed at mathematically demonstrating the effect of viscosity.
The molecular clutches bind between the actin filaments and the substrate with a reference association rate k 0 on , or equivalently at a timescale t 0 on ; 1=k 0 on .At any instant, the number of the associated clutches is n c ; PN 0 c .The substrate force builds up at an FA site with the actin filaments moving toward the cell nucleus with the retrograde speed ðV r ; _ x r Þ.The total restoring force of the bounded clutches is balanced by the substrate deformation.By first considering an elastic substrate, we have K e x sub = n c K c (x r − x sub ).Since n c ≫ 1, it is easy to show that the clutch deformation (x r − x sub ) is insignificant compared with the substrate deformation x sub , or x r ≈ x sub [45].Therefore, we have the force-velocity relation in an alternative form, On an elastic substrate, equation (A 3) only contains one unknown x sub and can be solved with the initial condition x sub | t=0 = 0, The substrate force is thus given by The timescale t l ¼ N 0 m fm V0Ke sets the substrate deformation rate and can approximate the lifetime of a binding-unbinding cycle in the FA dynamics [28].By comparing the substrate deformation timescale τ l with the binding timescale t 0 on , we defined a threshold stiffness K 0 ¼ N m f m k 0 on =V 0 (table 2).On soft substrates with K e /K 0 < 1, the clutch tension builds up slowly, and a long lifetime ðt l .t 0 on Þ allows for the association of a large number of clutches, i.e. n c N 0 c .On stiff substrates, the traction force increases rapidly, and a short lifetime ðt l , t 0 on Þ limits the clutch association.Consequently, the force of a single engaged clutch increases significantly (f > f cr in equation (2.3)), triggering adhesion reinforcement.
By describing a viscoelastic substrate by the SLS model, the substrate force-displacement relation is given by Since equation (A 3) still holds on viscoelastic substrates as long as the instantaneous stiffness is below the threshold, i.e.K e + K a < K 0 , we substitute equation (A 3) into equation (A 6) to eliminate F sub , which gives where τ r ≡ γ/K a defines the relaxation timescale.With the initial conditions x sub = 0 and F sub = 0 at t = 0, the solution of the substrate force is obtained as where τ 1 and τ 2 represent two characteristic times with τ 1 ≡ τ l + τ r (1 + (K a /K e )), and τ 2 ≈ τ r τ l /τ 1 .For soft substrates with a low equilibrium stiffness K e ≪ K a (i.e.K a /K e ≫ 1), we have τ 2 /τ l = 1/((1 + K a /K e ) + τ l /τ r ) ≪ 1 and τ 2 /τ 1 = 1/(2(1 + K a / K e ) + τ l /τ r + τ r /τ l (1 + K a /K e ) 2 ) ≪ 1.We can thus simplify the equation as The solutions have good agreement with the numerical results prior to the clutch rupturing at t ≈ τ l (electronic supplementary material, figure S15B).Equation (A 9) can be simplified if we focus on two extreme cases.For a substrate with a short relaxation time ðt r t 0 on Þ, the long duration of the binding-unbinding cycles allows us to neglect the influence of short characteristic time τ 2 , which gives where t 2 % N 0 m f m =ðK a þ K e Þ, same as the characteristic time on an elastic substrate (equation (A5)) with the stiffness K a + K e .
With these analytical solutions, we can calculate the average substrate force of binding/unbinding cycle by F sub ¼ Ð t l 0 F sub dt=t l .When t r t 0 on , the clutch system preserves the lifetime on an elastic substrate, i.e. t l % N 0 m f m =V 0 K e .By contrast, a long relaxation time means that the substrate deforms like an elastic substrate with a stiffness K a + K e .As a result, the shortened lifetime is comparable to the characteristic timescale τ 2 of the substrate deformation when τ r > τ l .Thus, the average substrate force in two limits is given by , i ft r !t l : Although the average substrate traction force is difficult to solve for when t 0 on , t r t l , our numerical results demonstrate a fast decrease of the clutch lifetime when t r .t 0 on , leading to a quick drop of the average substrate force (electronic supplementary material, figure S15C).Taken together, the average traction force must vary in a biphasic fashion with increasing viscosity and achieve a maximum when t r % t 0 on (electronic supplementary material, figure S15D).The above derivations approximate the dynamics of the traction force at a single FA site with N c ≈ N m .We next briefly explain the impact of viscosity on the net traction force, which is the traction difference between the cell front and rear.On soft substrates with K a + K e < K 0 , the traction force decreases with increasing viscosity beyond g .t 0 on K a due to the reduced FA lifetime (electronic supplementary material, figure S15C).Consequently, the net traction force must also vary biphasically with increasing viscosity.Owing to the regulation of Rho GTPase molecules, the front of the cell has a higher number of clutches than motors (N i c .N i m ), reducing the tension (f i c ) shared by each bond and thus prolonging the clutch lifetime (electronic supplementary material, figure S15E).The extended FA lifetime preserves the high traction force at the front of the cell even when t r .t 0 on (electronic supplementary material, figure S15F).This explains why the viscosity that maximizes the net traction force is slightly higher than the viscosity that produces the maximum traction force at the rear FA sites.

A.3. Estimation of individual clutch force under adhesion reinforcement
On very stiff substrates, the rapid building of tension f i c on an individual clutches contributes to a large average bond force far beyond the threshold force (i.e.hf a i t l ) f cr ), giving rise to very strong adhesion reinforcement effect.The association rate can thus be approximated by k i on % k 0 on exp ½zðhf a i t l À f cr Þ.For an engaged clutch with well-developed tension f i c , the clutch disassociation rate is mainly dominated by effect of the slip bond, i.e. k i off % k r0 expð f i c =f r0 Þ.For the saturation of the bounded clutches on very stiff substrates, we must have dP i /dt = 0 and k i on ¼ ð1 À P i Þ=P i Þk i off % ck i off (equation (2.4)), where c denotes a constant.Therefore, we know that the bond force sustained by an individual clutch can be related to the adhesion reinforcement coefficient by f i c =f r0 / zðhf a i t l À f cr Þ.Typically, the maximum of hf a i t l in our numerical analyses is around 10-12 pN.We can thus further obtain that f i c =f r0 / 8z pN.

Figure 1 .
Figure 1.Multiscale whole-cell theory schematics.(a) Intracellular organization involved in mesenchymal migration.(b) Close-up view of the cytoskeleton-ECM linkage via integrins and the adaptor proteins (e.g.talin).(c) The parameters and variables of the motor-clutch model with adhesion reinforcement at each FA site.(d ) The vertex-based model that couples the cell-matrix interactions at the subcellular scale to the chemomechanical dynamics at the cellular scale.The Rac1 concentration R and RhoA concentration ρ can each be in a membrane-bound active (red dots) or inactive state (orange squares), or dissolved in the cytoplasm (grey triangles).The conversion rates between these states areA i G , I i G , M À G , M þ G (G ≡ {R, ρ}); D is the diffusion constant of the membranebound species (equation(2.15)).The angle between the two membrane sections at a vertex is denoted by θ i .(e) The principal dimensions l max , l med , l min of a deformed nucleus.

Figure 2 .
Figure2.Cell spreading area and migration speed on elastic substrates.(a) The dimensionless mean cell area A=L 2 (full line) as a function of the dimensionless substrate elastic stiffness K e /K 0 .For comparison, the dots and error bars are the experimental data from[6].1(b) The average dimensionless substrate traction force versus K e /K 0 at vertex 0 (cell front) and vertex 8 (cell rear).The red dash-dotted line gives the upper bound of the substrate traction force on very stiff substrates.The calculation of the averaged substrate traction forces hF i sub i t total from the full focal adhesion cycles is illustrated in electronic supplementary material, figure S11.(c) The dimensionless mean and standard deviations of the net traction force versus K e /K 0 .The calculation of the time-averaged net traction force is illustrated in electronic supplementary material, figure S12.(d ) Cell-spreading induced nucleus viscous drag increase f shape h cp =h 0 cp versus K e /K 0 , which is calculated by using equations (2.9), (2.11), (2.12) for the mean spreading area values in a. (e) The simulated dimensionless mean migration speed V=V 0 versus K e /K 0 (full curve) and the experimental data for the U373-MG human glioma cells (dots and error bars)[6].( f ) The influence of the shape factor f shape of the deformed nucleus and the area-dependent cytoplasm viscosity η cp on the migration speed (equations (2.11), (2.12)).The mean values or the standard deviations are calculated over n = 10 simulations at each data point in a, c-f.

Figure 4 .
Figure 4. Simulations reproduce cell durotaxis on elastic gradient substrates.(a) Trajectories of 30 simulated cells over greater than 7.2 h in real units on an elastic substrate with a linear stiffness variation between K e = 0.5-22 pN nm −1 in the y-direction.(b) The percentage of the simulated cells out of n = 83 simulations on the bottom, middle and top third of the gradient substrate as a function of time, i.e. at t = 0 h (initial condition), t = 7.2 h, t > 12 h, which are compared with the distributions of the MDA-MB-231 cells (grey dashed lines) that are seeded on polyacrylamide hydrogels with a stiffness gradient and observed for 72 h[16].

Figure 5 .
Figure 5. Directed migration is regulated by viscosity gradient on soft substrates.Representative trajectories of 20 cells and angular displacements of approximately 50 cells on a viscoelastic substrate (K a = 0.1 pN nm −1 , K e = 2.9 pN nm −1 ) with (a), (b) a linearly varying viscosity γ = 0.0-1 pN s nm −1 , (c), (d ) γ = 1.0-20 pN s nm −1 , and (e), ( f ) a uniform viscosity γ = 10 pN s nm −1 .(g)The simulated versus analytical traction forces (averaged over τ l ) as a function of the ratio of the material relaxation timescale τ r ≡ γ/K a to the clutch binding timescale t 0 on ; 1=k 0 on (table2).The traction forces at the two opposite ends of a cell along the direction of the motion are labelled by F R am or F L am , respectively.

Table 1 .
Residual substrate stiffness K e , additional stiffness K a and viscosity γ ranges in simulations.
.org/journal/rsif J. R. Soc.Interface 20: 20230317 second term in the square bracket negligible, i.e.F sub % N 0 m f m 1 À exp À t t 2, when t r !t l , ðA 1 À on : ðA 10Þ By contrast, a very viscous substrate corresponds to a long relaxation timescale τ r ≥ τ l , making the contribution of the royalsocietypublishing